#rm(list=ls())
library(boot);library(popbio);library(rjags);library(coda)
root = "/Users/Tom/Documents/YNP_Bison_Model/"

path = function (stem,r=root){
return(paste(r,stem, sep=""))	
}

load(path("A_final_model_runs_for_revised_manuscript/Frequency Dependent/One_beta_final.Rdata"))

setwd(path("A_final_model_runs_for_revised_manuscript/Analyses/"))

set.seed(1)

b=One_beta.coda
c=as.data.frame(rbind(b[[1]],b[[2]],b[[3]]))

A=matrix(0,8,8)

I=diag(8)
t=1

N.eq = 3000  #assumed size of healthy population including males--reults are very insensitive to this number. 





R0_f=function(rep=3000,N.eq=3000,v_switch=1,psi_switch=1,t=1){
p=numeric(3)
phi=numeric(1)
f.con=numeric(rep)
f.neg=numeric(rep)
f.pos=numeric(rep)
 b=numeric(3)
R0=real(rep)

for(j in 1:rep){
	phi=numeric(1)
	b=numeric(4)

	i=round(runif(n=1, min=1, max= nrow(c)))
	
	#get parameter estimates from MCMC chain

	b[1]=c$"b[1]"[i]
	b[2]=c$"b[2]"[i]
	b[3]=c$"b[3]"[i]
	#b[4]=c$"b[4]"[i]
	beta=c$beta[i]
	m=c$m[i]
	p[1]=c$"p[1]"[i]
	p[2]=c$"p[2]"[i]
	p[3]=c$"p[3]"[i]
	psi=c$psi[i]*psi_switch
	v=c$v[i]*v_switch
	f.pos=c$"b[2]"[i]
	f.neg = c$"b[1]"[i]
	f.con = c$"b[3]"[i]
	
	
	#cacluate stable stage composition of healthy population
	M=matrix(0,nrow=4,ncol=4)
	M[1,3]=f.neg
	M[2,1]=p[1]*m
	M[3,2]=p[2]
	M[3,3]=p[2]
	M[4,1]=p[1]*(1-m)
	M[4,4]=p[3]
	e=eigen.analysis(M)
	n=e$stable.stage*N.eq
	#cacluate probabilty that a susecptible becomes infected when there is a single infected in stage 6. Note that sum(n[1:3]) is the total number of females and young in the population.  Note t is set as a constant = 1.
	phi[t]=(1-exp(-c$beta[i]*p[2]/sum(n[1:3])))   
	
A[1,3]<-p[2]*f.neg[t]*(1-phi[t])
A[1,6]<-p[2]*f.con[t]*(1-v)*(1-phi[t])
A[1,7]<-p[2]*(f.pos[t] *(1-psi)*(1-phi[t])  +  f.pos[t] * psi*(1-v)*(1-phi[t]))  #changes


#row 2
A[2,1]<-p[1]*(1-m)*(1-phi[t])
#row 3
A[3,2]<-p[2]*(1-phi[t])
A[3,3]<-p[2]*(1-phi[t])

#row4 
A[4,3]<-p[2]*f.neg[t]*phi[t]
A[4,6]<-p[2]*f.con[t]*(v+phi[t]-v*phi[t])
A[4,7]<-p[2]*(f.pos[t]*psi*(v+phi[t]-v*phi[t]) + f.pos[t]*(1-psi)*phi[t])

#row5
A[5,1]<-p[1]*phi[t]*(1-m)
A[5,4]<-p[1]*(1-m)

#row6
A[6,2]<-p[2]*phi[t]
A[6,3]<-p[2]*phi[t]
A[6,5]<-p[2]
A[6,7] <- 0

#row 7
A[7,6]<-p[2]
A[7,7]<-p[2]


#row 8
A[8,1]<-p[1]*m
A[8,4]<-p[1]*m
A[8,8]<-p[3]

 	T=A
	T[1,]=0
	T[4,]=0
	N=solve(I-T)
	N
	s=p
	F=matrix(0,8,8)
		
	
	F[6,6] = phi[t]*p[2]*(n[3] + n[2] )#The number of new infected individulas in stage j (6) created by an indvidual in stage i (6)
	 #The number of new infected juveniles contributed by an individual in stage 6.  Multiply by 1-m because only considering offspring that can survive to become infectious
	F[4,6]<-f.con[t]*(v+phi[t]-v*phi[t])*(1-m)
	#The number of new infected yearlings contributed by an individual in stage 6
	F[5,6] = phi[t]*p[1] * n[1] 
	#The number of new infected indviduals created by an indvidual in stage 7

	F[6,7] <- psi*phi[t]*p[2]*(n[3]+n[2])
	#The number of new infected juveniles contributed by an individual in stage 7
	F[4,7] = f.pos*(psi*(v+phi[t]-v*phi[t]) + (1-psi) * phi[t])*(1-m)
	#The number of new infected yearlings created by an individual  in stage 7
	F[5,7] =  psi*p[1]*phi[t]*n[1] 
	
	
	
	
	R=F %*% N
	ev=eigen(R)
	
	R0[j]=ev$values[1]

} #end of loop
	return(as.real(R0))
}#end of function

rep=50000
f=R0_f(rep=rep)
no_v=R0_f(rep=rep, v_switch=0)
no_psi=R0_f(rep=rep, psi_switch=0)

R0=list(f=f,no_v=no_v,no_psi=no_psi)

#plot(density(R0),xlim=c(0,20))
par(mfrow=c(1,1))
hist(f,breaks=200,freq=FALSE,xlab=expression(paste("Basic reproductive ratio, ", R[0])),xlim=c(0,10))

summary(f)
quantile(f, c(.025, .975), na.rm=TRUE)
1-ecdf(f)(5)




save(R0, file= "R0.data.Rdata")
#